//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : RrbSolver.h
//----------------------------------------------------------------------------


#ifndef RRBAKASHSOLVER_H
#define RRBAKASHSOLVER_H
#include <vector>
#include <cmath>
#include "Array2D.h"
#include "defines.h"
#include <omp.h>

//#define prec
//----------------------------------------------------------------------------
// class RrbSolver
//----------------------------------------------------------------------------
/// recursive red black solver by Auke van der Ploeg, Coupled with BiCGSTAB by Akash
//----------------------------------------------------------------------------
class RrbAkashSolver
{

public:

    // Constructor: Pass the 9 Diagonals required.
    RrbAkashSolver(const Array2D<REAL> &pzzcc,
                   const Array2D<REAL> &pzzcs,
                   const Array2D<REAL> &pzzcw,
                   const Array2D<REAL> &pffcc,
                   const Array2D<REAL> &pffcs,
                   const Array2D<REAL> &pffcw,
                   const Array2D<REAL> &pzfcc,
                   const Array2D<REAL> &pzfcs,
                   const Array2D<REAL> &pzfcw,
                   const Array2D<REAL> &pfzcc,
                   const Array2D<REAL> &pfzcs,
                   const Array2D<REAL> &pfzcw,
                   int Nx1,
                   int Nx2,
                   REAL tolerance);

    virtual ~RrbAkashSolver();

    void solve(Array2D<REAL> *pB1,
               Array2D<REAL> *pB2,
               Array2D<REAL> *pX1,
               Array2D<REAL> *pX2);
    void GMRES(Array2D<REAL> *pB1,
               Array2D<REAL> *pB2,
               Array2D<REAL> *pX1,
               Array2D<REAL> *pX2);


private:

    int m_Nx1;
    int m_Nx2;
    int COLUMNS;
    int ROWS;
    REAL m_tolerance;
    /// @cond nevercall
    RrbAkashSolver(const RrbAkashSolver& m);
    RrbAkashSolver & operator =(const RrbAkashSolver& m);
    /// @endcond


    void preconditioner(const std::vector<Array2D<REAL> > &mat_st,std::vector<Array2D<REAL> > *mat_pc);
    void solver(const std::vector<Array2D<REAL>> &mat, Array2D<REAL> *p_cse );
    void solver2(const std::vector<Array2D<REAL>> &mat, Array2D<REAL> *p_cse);

    void storematrices();
    void matvec_rrb(const Array2D<REAL> &X1, const Array2D<REAL> &X2 , Array2D<REAL> *Ax1, Array2D<REAL> *Ax2);
    void prec_solve(Array2D<REAL> &Xz_hat, Array2D<REAL> &Xf_hat);
    void prec_solve_2(Array2D<REAL> &Xz_hat, Array2D<REAL> &Xf_hat);
    void prec_solve_3(Array2D<REAL> &Xz_hat, Array2D<REAL> &Xf_hat);


    inline void GeneratePlaneRotation (const REAL &dx, const REAL &dy, REAL *c1, REAL *s2);
    inline void ApplyPlaneRotation (REAL *dx, REAL *dy,const REAL &cs,const REAL &sn);
    inline void Update (Array2D<double> *X1, Array2D<double> *X2,
                        int k, const Array2D<REAL> &h,
                        const std::vector<REAL> &s,
                        const std::vector<Array2D<REAL> > &Vz,
                        const std::vector<Array2D<REAL> > &Vf);
    void ComputeResidual(const Array2D<REAL> &Xz,const Array2D<REAL> &Xf,
                                         const Array2D<REAL> &Bz,const Array2D<REAL> &Bf,
                                         Array2D<REAL> *R1 , Array2D<REAL> *R2);


/*
    REAL dotproduct(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2);
    REAL norm(const Array2D<REAL> &Vec1,const Array2D<REAL> &Vec2);
    REAL norm(const Array2D<REAL> &Vec1);


    Array2D<REAL> MatAdd(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2);
    Array2D<REAL> MatSub(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2);
    Array2D<REAL> scalarmult(const REAL &mult, const Array2D<REAL> &Vec);
*/



    Array2D<REAL> m_pzzcc; /// pointer to MatzC - ZETA
    Array2D<REAL> m_pzzcs; /// pointer to MatzS   ZETA
    Array2D<REAL> m_pzzcw; /// pointer to MatzW   ZETA

    Array2D<REAL> m_pffcc; /// pointer to MatfC   PHI
    Array2D<REAL> m_pffcs; /// pointer to MatfS   PHI
    Array2D<REAL> m_pffcw; /// pointer to MatfW   PHI

    Array2D<REAL> m_pzfcc; ///
    Array2D<REAL> m_pzfcs; ///
    Array2D<REAL> m_pzfcw; ///

    Array2D<REAL> m_pfzcc; ///
    Array2D<REAL> m_pfzcs; ///
    Array2D<REAL> m_pfzcw; ///

    Array2D<REAL> m_czzne; /// K1
    Array2D<REAL> m_czznw; /// K2
    Array2D<REAL> m_czzse; /// Az
    Array2D<REAL> m_czzsw; /// K4

    Array2D<REAL> m_cffne; /// K1
    Array2D<REAL> m_cffnw; /// K2
    Array2D<REAL> m_cffse; /// Az
    Array2D<REAL> m_cffsw; /// K4

    Array2D<REAL> m_czzcc; /// copyMatC
    Array2D<REAL> m_czzss; /// copyMatS
    Array2D<REAL> m_czzww; /// copyMatW
    Array2D<REAL> m_czzee; /// copyMatE
    Array2D<REAL> m_czznn; /// copyMatN

    Array2D<REAL> m_cffcc; /// copyMatC
    Array2D<REAL> m_cffss; /// copyMatS
    Array2D<REAL> m_cffww; /// copyMatW
    Array2D<REAL> m_cffee; /// copyMatE
    Array2D<REAL> m_cffnn; /// copyMatN

    // Vectors storing the pentadiagonals
    std::vector<Array2D<REAL> > m_ZetaZeta;
    std::vector<Array2D<REAL> > m_PhiPhi;
    std::vector<Array2D<REAL> > m_ZetaPhi;
    std::vector<Array2D<REAL> > m_PhiZeta;

    // Vectors storing the Preconditioners
    std::vector<Array2D<REAL> > m_pZetaZeta;
    std::vector<Array2D<REAL> > m_pPhiPhi;



    int m_st;
    int m_tw;
    int m_niv;

    REAL m_rhobew;
    REAL m_bet;
    REAL m_alf;
    bool m_done;
    bool m_done1;
    bool m_done2;
};
/*
inline Solver* NewRrbSolver(Array2D<REAL> *pmatzC, ///< central matrix
                            Array2D<REAL> *pmatzS, ///< south matrix
                            Array2D<REAL> *pmatzW, ///< west matrix
                            Array2D<REAL> *pmatfC, ///< central matrix
                            Array2D<REAL> *pmatfS, ///< south matrix
                            Array2D<REAL> *pmatfW, ///< west matrix
                            Array2D<REAL> *pmatzfC, ///< central matrix
                            Array2D<REAL> *pmatzfS, ///< south matrix
                            Array2D<REAL> *pmatzfW, ///< west matrix
                            Array2D<REAL> *pmatfzC, ///< central matrix
                            Array2D<REAL> *pmatfzS, ///< south matrix
                            Array2D<REAL> *pmatfzW, ///< west matrix
                            int N, ///< width
                            int M, ///< height
                            REAL tolerance) ///< solver tolerance
{
    return new RrbAkashSolver(pmatzC, pmatzS, pmatzW,pmatfC, pmatfS, pmatfW,pmatzfC, pmatzfS, pmatzfW,pmatfzC, pmatfzS, pmatfzW, N, M, tolerance);
}*/

#endif // RRBSOLVER_H
